function [CE,P_star] = Get_CE_LSF_optim_fun(P_choices,MF_fee_rate,u,T,g_I,R_S_GHQ,weight,r,A_grid,I_0)

%% Set model specifications & parameters
P_num = length(P_choices);


%% State space
A_num = length(A_grid); A_grid_max = max(A_grid);   % grid of potential QRP account values A_t.


%% Matrizes to store results from dynamic optimization
EU_at = zeros(A_num,T);           % to store expected utility of being in current state (and time)


%% Terminal utility
EU_at(:,T) = u(max(A_grid,0.01));


%% Recursion
EU_0 = zeros(P_num,1);
for p = 1:P_num
    P = P_choices(p);

    R_A = max( P * R_S_GHQ + (1-P) * (exp(r)-1) , -1 );             % GHQ nodes for return credited to QRP account.
    
    for t = (T-1):-1:1
        EU_tplus1 = EU_at(:,t+1);     % value function at end of period.
        I_t = I_0 * (1+g_I)^t;  % PH's time-t contribution to QRP (no contribution at retirement time T).
        for a = 1:A_num
            A_t = A_grid(a);
            A_tplus1 = min( (A_t + I_t)*(1+R_A)*exp(-MF_fee_rate) , A_grid_max );         % GHQ nodes for QRP account value at year-end.
            U_tplus1 = interp1(A_grid,EU_tplus1,A_tplus1);              % GHQ nodes for utility at year-end.
            EU_at(a,t) = weight' * U_tplus1;                               % expected utility of account values at end of period.
        end
    end

    % Time 0
    EU_tplus1 = EU_at(:,1);     % value function at end of period.
    I_t = I_0;                  % PH's time-t contribution to QRP (no contribution at retirement time T).
    A_t = 0;
    A_tplus1 = min( (A_t + I_t)*(1+R_A)*exp(-MF_fee_rate) , A_grid_max );         % GHQ nodes for QRP account value at year-end.
    U_tplus1 = interp1(A_grid,EU_tplus1,A_tplus1);              % GHQ nodes for utility at year-end.
    EU_0(p) = weight' * U_tplus1;                               % expected utility of account values at end of period.

end
[EU_star_0,p] = max(EU_0);
P_star = P_choices(p);


%% Convert utility to Certainty Equivalent (CE)
options = optimset('Display','none');
CE = fsolve(@(ce)  1e+3 * (u(ce)/EU_star_0-1), 200,options);    % certainty equivalent to potential values of A_T.

return
